﻿using LatoolNet;

namespace NumericalSolution {
  class Diffusion1DImplicit : Diffusion1D {
    public override double[] Next() {

      base.Age();

      int rownum = p_numberOfNodes;
      int colnum = p_numberOfNodes;
      int bandWidth = 3;
      Matrix mat = new Matrix(rownum, colnum, bandWidth);
      Matrix vec = new Matrix(colnum, 1);

      double s = (p_alpha * p_deltaTime) / (p_deltaX * p_deltaX);

      mat[0, 0] = 1;
      vec[0, 0] = 0.0;
      for (int i = 1; i < p_numberOfNodes - 1; i++) {
        mat[i, i - 1] = -s;
        mat[i, i] = (1 + 2 * s);
        mat[i, i + 1] = -s;
        vec[i, 0] = p_previous_u[i];
      }
      mat[p_numberOfNodes - 1, p_numberOfNodes - 1] = 1;
      vec[p_numberOfNodes - 1, 0] = 0.0;

      //LUFactorization.Solve(mat, vec);
      LinearEquation.Solve(mat, vec);

      for (int i = 0; i < p_numberOfNodes; i++) {
        p_current_u[i] = vec[i, 0];
      }

      mat.Dispose();
      vec.Dispose();
      return p_current_u;
    }
  }
}
